﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.NumericalIntegration
{
    /// <summary>
    /// This class provides the integration of a complex integral by using the Exponential Transform method.
    /// </summary>
    [Serializable]
    public class ComplexExponentialTransformIntegrator 
    {
        /// <summary>
        /// This field holds the integral to integrate.
        /// </summary>
        private IComplexFinitelyIntegral integral;

        /// <summary>
        /// Initializes a new instance of the <see cref="ComplexExponentialTransformIntegrator"/> class.
        /// </summary>
        /// <param name="integral">The integral for the numerical integration.</param>
        public ComplexExponentialTransformIntegrator(ComplexFinitelyIntegral integral)
        {
            this.integral = integral;
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="ComplexExponentialTransformIntegrator"/> class.
        /// </summary>
        /// <param name="integral">The integral for the numerical integration.</param>
        public ComplexExponentialTransformIntegrator(IComplexFinitelyIntegral integral)
        {
            this.integral = integral;
        }

        /// <summary>
        /// Gets or sets the integral to integrate.
        /// </summary>
        /// <value>The integral to integrate.</value>
        public IComplexFinitelyIntegral Integral
        {
            get { return this.integral; }
            set { this.integral = value; }
        }

        /// <summary>
        /// Integrates the integral in a numerical way.
        /// </summary>
        /// <param name="iterations">The maximum numbers of iterations.</param>
        /// <returns>The area of the specified integral.</returns>
        public ComplexNumber Integrate(int iterations)
        {
            iterations = iterations / 2;

            const double pi4 = Math.PI / 4;
            ComplexNumber h = new ComplexNumber(5.0 / iterations);
            ComplexNumber ss = 0.5 * this.integral.Function.SolveAt((this.Integral.A + this.Integral.B) / 2.0);

            for (int k = -iterations; k < 0; k++)
            {
                ComplexNumber z = h * new ComplexNumber(k);
                ComplexNumber exz = ExtendedMath.Exp(z);
                ComplexNumber hcos = exz + 1.0 / exz;
                ComplexNumber hsin = exz - 1.0 / exz;
                ComplexNumber s = ExtendedMath.Exp(pi4 * hsin);
                ComplexNumber w = s + 1.0 / s;
                ComplexNumber dxdz = hcos / (w * w);
                ComplexNumber x1 = (this.Integral.B * s + this.Integral.A / s) / w;
                ComplexNumber x2 = (this.Integral.A * s + this.Integral.B / s) / w;

                if ((x1 != this.Integral.A) && (x1 != this.Integral.B))
                {
                    ss += dxdz * this.integral.Function.SolveAt(x1);
                }

                if ((x2 != this.Integral.A) && (x2 != this.Integral.B))
                {
                    ss += dxdz * this.integral.Function.SolveAt(x2);
                }
            }

            return 2.0 * (this.Integral.B - this.Integral.A) * pi4 * h * ss;
        }
    }
}
